clc; clear; close all;
% 固定系统参数 mu=3.5、初值 (0.1,0,−1) 的条件下，NMSE 随耦合强度 k 变化的情况


% 参数设置
mu = 3.5;
x0_init = 0.1; 
y0_init = 0;
z0_init = -1;        % 固定 (x0, y0, z0)

nk = 1000;            % k 的采样数
k_vals = linspace(0, 0.2, nk);  % 扫描区间 [0, 0.2]

% 迭代步数设置
N_transient = 500;   % 瞬态迭代步数
N_nmse      = 500;   % 用于计算 NMSE 的采样步数

NMSE_array = zeros(1, nk);  % 用于存储每个 k 的 NMSE

% 扫描 k 计算 NMSE
for i = 1:nk
    k = k_vals(i);

    % 1) 先做瞬态迭代，去除初始影响
    state = [x0_init, y0_init, z0_init];
    for n = 1:N_transient
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
    end
    
    % 2) 从瞬态后的状态出发，正式计算 NMSE
    sum_err = 0; 
    s = state;
    for n = 1:N_nmse
        [dx, dy, dz] = mclm(s, mu, k);
        s = [dx, dy, dz];
        
        denom = sqrt(dx^2 + dy^2);
        if denom < 1e-12
            % 分母过小时可跳过或设为阈值
            continue;
        end
        sum_err = sum_err + abs(dx - dy)/denom;
    end
    
    % 取平均，得到 NMSE
    NMSE_array(i) = sum_err / N_nmse;
end

% 绘 NMSE 曲线
figure;
plot(k_vals, NMSE_array, 'g-', 'LineWidth',1);
xlabel('k');
ylabel('NMSE');
title('NMSE 与 k 的关系');
grid on;
